Intrinsic Defects and Dopability of Zinc Phosphide 
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Zinc Phosphide {ZnzP2) could be the basis for cheap and highly efficient solar cells. Its use in 
this regard is limited by the difficulty in n-type doping the material. In an effort to understand 
the mechanism behind this, the energetics and electronic structure of intrinsic point defects in zinc 
phosphide are studied using generalized Kohn-Sham theory and utilizing the Heyd, Scuseria, and 
Ernzerhof (HSE) hybrid functional for exchange and correlation. Novel 'perturbation extrapolation' 
is utilized to extend the use of the computationally expensive HSE functional to this large-scale 
defect system. According to calculations, the formation energy of charged phosphorus interstitial 
defects are very low in n-type Zn3P2 and act as 'electron sinks', nullifying the desired doping 
and lowering the fermi- level back towards the p-type regime. This is consistent with experimental 
observations of both the tendency of conductivity to rise with phosphorus partial pressure, and with 
current partial successes in n-type doping in very zinc-rich growth conditions. 
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I. INTRODUCTION 

Zinc Phosphide has great potential as a photovoltaic 
material; it absorbs strongly in the visible spectrum 
(> 10^ — 10^c?n~^), has long minority carrier diffusion 
lengths (5-10 /iTTi), and a direct, almost ideal bandgap 
(~1.5eV) % Both zinc and phosphorus are Earth- 
abundant elements, greatly aiding their widespread use. 
However, no practical means of creating Zn^P2 crystals 
with n-type doping has been found. This has prevented 
the typical p-n homo junction solar cells and current 
implementations instead rely upon metal-semiconductor 
junctions or p-n semiconductor hetero junctions^. The 
best results are currently with p-Zn3P2/Mg Schottky 
diode where the maximal open-circuit voltage (~0.5eV) 
is the limiting factor on the efficiency of these devices 
(currently 6%)"^. 

Differences in growth environment has been shown to 
affect the electrical properties significantly^. Specifically, 
resistivity ranging over 100 to lO^l^-cm have been mea- 
sured in single and polycrystalline samples. Such a large 
range of values may point to an intrinsic defect mech- 
anism dominating carrier concentrations. Despite many 
years of research, the exact nature of the various intrinsic 
defects as well as their small-scale structure and proper- 
ties are only partially known. So, as a step towards un- 
derstanding the n-type doping difficulties, it is natural to 
first explore the role of intrinsic defects in this system. 

Density Functional Theory (DFT) is the method of 
choice for studying the electronic structure of defects in 
semiconductors. It is the only ah initio approach cur- 
rently tractable for calculation of the energetics of cells 
of sufficient size to explore isolated defects (at least on 
the order of lOO's of atoms). For computational sim- 
plicity, the most common exchange-correlation function- 
als utilized are the local-density approximation (LDA) 
and generalized-gradient approximation (GGA). How- 
ever, while the defect-cell total energies are typically 
well- described by LDA or GGA, the well-known bandgap 
problem of these methods poses problems in the de- 



termining the precise electronic structure of the defects 
themselves. Subsequently, this can lead to different pre- 
dictions of defect stability amongst calculations even on 
identical systems. For example, recent work of various 
groups on ZnO has lead to predictions that the oxygen va- 
cancy acts as both a shallow and deep defect^ -. Recently, 
more accurate techniques such as hybrid functionala^ and 
GW excited-state calculations^ yield greatly improved 
bandgap prediction. However, these techniques are still 
too expensive for the large scales required of defect sys- 
tems. 

Regardless of the exact ah initio procedure chosen, 
there is a fundamental problem in studying defects with 
a necessarily small amount of atoms (even the most ef- 
ficient approach of LDA-based DFT is limited to on the 
order of hundreds of atoms); that is, we are imposing a 
degenerate doping condition, which is seldom the regime 
of interest (typical defect or dopant concentrations rarely 
exceed parts-per-million). As such, care must be taken 
to correct for interactions that are artifacts of this incor- 
rectly high defect concentration such as image charges 
and spurious hybridization. In the following section we 
detail the methods that allow us to compute defect levels 
and energetics separate from these effects. 

ZnsP2 exists in two phases, tetragonal and cubic. 
We have investigated the intrinsic point defects of the 
tetragonal phase, which is of primary concern because 
it is the room-temperature phase. These defects include 
zinc vacancy (Vzn), phosphorus vacancy (Vp), zinc 
interstitial (Zn^), and phosporus interstitial {Pi). The 
antisite defects were also studied, but their formation 
energies are such that they would be unstable to dis- 
sociating into a vacancy-interstitial pair, thus they are 
excluded from the following discussion. What we find 
is that the formation energy of the charged acceptor 
defects, especially P^, is small enough at even moderate 
n-type conditions to form enough compensating defects 
as to completely neutralize the desired extrinsic doping. 



II. METHODOLOGY 

Of central concern is the stability of the intrinsic de- 
fects in ZnsP2 and their effect on the fermi level. The 
stability of a defect is largely determined by its formation 
energy; in the supercell formalism put forth by Lany and 
Zungei^ the formation energy of a defect with charge q 
is given by a sum of three terms: 



AHf = [Ed - Eh] + q{Ev + AEf) 






Uai/J^a ^^l^o 



(1) 

Eh and Ed are the total energy of the perfect host and 
host+defect supercells respectively. The first term on the 
right is the difference in bond energy brought about by 
the defect. 

The second term in ^ represents the energy cost of 
exchanging electrons with the 'electron reservoir'. Ey 
is the reference energy of the reservoir and is the price 
we pay for removing an electron from the top of the va- 
lence band (ie. the energy of a hole at the Valence Band 
Maximum)?^. Consequently, equation Q describes the 
energetics of forming a defect while conserving charge^. 
The calculated total energy (Eh) of the system follows 
Janak's Theoremffl 



drij 



= ei 



(2) 



where Ui is the occupation of the highest occupied state i 
with eigenvalue e^ . For an infinite system e^ is identical to 
Ey. Thus, Ey can be calculated as the energy difference 
between a host and a host+hole cell in the limit that the 
number of electrons (N) tends to infinity: 



lim [Eh{N) - Eh{N ■ 



1)] = Ev 



(3) 



As a practical matter, a good approximation can be at- 
tained for relatively small systems - in the present work 
the difference between a neutral perfect supercell and a 
supercell and a hole {Eh{0) - Eh{-\-)) is used. We now 
have a good approximation for the electron chemical po- 
tential at the VBM. Finally, AEf is the additional en- 
ergy of electrons in our system above the VBM and is 
the proxy for specifying the doping regime (p or n-type) 
of the bulk. 

The crystal growth environment affects the formation 
energy via the chemical potentials (/i^ = Ma + ^Ma) in 
equation ^ . These represent the energy cost of exchang- 
ing atoms with the chemical reservoir. By convention, 
the formation energies are defined in relation to the stan- 
dard states of the constituents of a system. Thus, the 
elemental chemical potentials are broken into two parts: 
/i^ is the chemical potential of the standard state of the 
element and A/Xq, is the chemical potential of the element 
in relation to the standard state. The growth conditions 



FIG. 1. Zn2,P2 allowed chemical potentials. The red and 
green lines are the stability limits for Zn^P2 and ZnP2 re- 
spectively. The dark purple region corresponds to the range 
of chemical potentials where we are assured of forming zinc 
phosphide crystals. 
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are reflected in A/iQ, (ie. a maximally rich growth en- 
vironment of a certain element would have A/iQ, = 0), 
it becomes more negative for lower concentrations of an 
element during crystal formation. The chemical poten- 
tials are added or subtracted from the formation energy 
according to the number of atoms of a certain species is 
deposited or withdrawn from the growth reservoir (no,); 
it is +1 if an atom is added (ie. vacancy defects), if we 
remove an atom it is -1 (ie. interstitial defects). 

Thus, the formation energy is determined as a func- 
tion of fermi- level {AEf) for a specific growth condition 
(dictated through A/Xq,). Assuming equilibrium growth 
conditions, the chemical potentials are restricted to val- 
ues that maintain a stable compound and don't permit 
competing phases to exist. For ZnsP2 we have the fol- 
lowing constraints: 

In order not to precipitate the elemental form of Zn or 
P, we must have: 



Mzn < Mzn; MP < Mp 



or, equivalent ly: 



A/izn < 0; Afip < 



(4) 



(5) 



where A/i^^=0 would define the maximumly-rich Zn 
growth condition. 



The stability condition for ZnsP2 is: 

3A/izn + 2A/ip < AHf{Zn^P2) 



(6) 



The only other competing phase is ZnP2 (phosphorus- 
rich phase) and avoiding its formation yields the following 
restriction: 



Ajizn + 2A/ip > AHf{ZnP2) 



(7) 



There is a the further issue of incomplete error cancella- 
tion in DFT when energy differences are taken between 



chemically dissimilar systems^ (such as between com- 
pounds and their elemental constituents). Here we cor- 
rect the chemical potentials according to [ 12 which re- 
sults in better agreement between predicted and exper- 
imental formation energies (1.04 eV non-corrected, 1.79 
eV corrected, 1.5 eV experimental). 

Combining all the conditions described by Eq 5-7, we 
can determine an allowed region of chemical potentials 
where we are assured of forming only Zn^P2 - see 
Figure [l] Later, when we discuss formation energies 
for a specific growth environment it will be this figure 
which defines the chemical potentials to use for the 
Zn-rich/P-poor versus the Zn-poor/P-rich regimes. 

All these calculations are performed within the Kohn- 
Sham framework (DFT) and utilize the projector aug- 
mented pseudopotentials as implemented in VASP^^Hi^. 
We use valence configurations of 3d4s and 3s3p for zinc 
and phosphorus respectively with the Perdew, Burke, 
and Ernzerhof (PBE) potential GGA as well as hybrid 
functionals (HSE) for exchange-correlation were used, 
the specifics of which will be discussed later. For super- 
cells of 2x2x2 unit cells (320 atoms) energies were calcu- 
lated with two k-points (0,0,0) and (1/2,1/2,1/2) and a 
plane-wave basis with a cutoff of 300eV. These settings 
result in a numerical error on the order of 0.1 eV in the 
defect formation energy which is on par with the other 
sources of error (e.g. uncertainty in the bandgap and 
finite-size corrections). All calculations were performed 
considering spin polarization. 

Although the use of a supercell geometry within a DFT 
framework is a common approach for defect calculations, 
care must be taken to avoid spurious finite-size effects 
as well as known deficiencies in DFT's ability to model 
excited state energies. 



1. Image Charge Correction 

The drawback of the use of a standard supercell geome- 
try is that defects are periodically and infinitely repeated 
spatially. The defect, instead of being surrounded by a 
large region of perfect bulk crystal as it would be under 
non-degenerate conditions, is now surrounded by mirror 
images of itself. This will result in somewhat frustrated 
ionic relaxation, though these elastic energy effects tend 
to be short range and is rarely a problem for even mod- 
estly sized cells (there is very little difference in relaxation 
energies for even a 2x2x2 supercell versus a unit cell of 
ZnsP2). However, when dealing with charged defects 
we form 'image charges' leading to spurious electrostatic 
interactions. These coulombic interactions between the 
defect and its mirror charges are long-ranged and signif- 
icant even for large cells. 

Corrections for this 'image charge' effect have been the 
subject of much research, though the most common ap- 
proach is based on the work of Makov and Payne^. They 
considered the charge density to be the contribution of 
the periodic charge of the underlying crystal structure 
and the charge density of the aperiodic defect (which is 
simply the electron density difference between the host 
and host+defect cells). The multipole correction to the 
formation energy is: 



Eic = 



2eL 



27rqQr 



■0{L-') 



(8) 



where aM is the supercell lattice-dependant Madelung 
constant, L is the length of the supercell, e is the static 
dielectric constant, and Qr is the second radial moment of 
the aperiodic charge density. The first two terms are the 
monopole and quadrupole corrections respectively; the 
quadrupole correction typically ~ 30% of the monopole 
terni^. 



Finite Size Effects 



With current limitations on computing power, DFT 
calculations are typically restricted to cells on the or- 
der of hundreds of atoms. Even the addition of a sin- 
gle point defect would thereby result in defect concen- 
trations on the order of tenths of a percent, normally 
describing degenerate conditions. Usual semiconductor 
defect concentrations are much more dilute (on the order 
of parts-per- million), and often result in far different ma- 
terial properties than within the degenerate regime. Due 
to the low dielectric strength (and hence poor screening) 
of ZnsP2 these effects are especially troublesome in this 
study. However, if we are careful about correcting our re- 
sults, we can still make accurate predictions about what 
the dilute environment should look like. The three main 
finite-size corrections are discussed below. 



2. Potential Alignment Correction 

In the case of charged defects with periodic bound- 
ary conditions there is a violation of charge neutrality, 
which causes the Coulomb potential to diverge^^. In 
momentum-space formalism, one usually sets the G=0 
term of the electrostatic and ionic potential {Vh{G=0) 
and Vi{G=0)) to zero. The Kohn-Sham eigenvalues are 
thus only defined with respect to the average electrostatic 
potential of the cell. For neutral systems this arbitrary 
offset still leads to a well-defined total energy since the 
electron-electron and ion-ion contributions exactly can- 
cel. In a charged system, ignoring the G=0 term can 
be viewed as equivalent to a uniform background charge 
(jellium) compensating for the net charge - though it is 
important to note that this only occurs for the potential. 
In a charged cell there is now an arbitrary offset to the 
total energy. The charged cell energies Ejj and Eh{-\-) 
in Eq. (fTj) have to be compensated for in order to treat 



FIG. 2. Potential alignment correction. Differences in elec- 
trostatic potential between perfect and defect cells far from 
defect area are used to correct charged cell energies. 
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them on an equal footing with the neutral cell^^. 
potential- alignment correction is: 



EpA = q{VL,, 



V^h) 



The 



(9) 



where the charged defect (Vg^) and host reference (V^) 
potentials are atomic sphere-averaged electrostatic po- 
tentials far from the defect site, q is the charge of the 
defect- see Figure [2] 



3. Band-filling Correction 

The artificially degenerate doping regime can cause de- 
fects to form bands rather than isolated states within 
the bandgap^^^^. For shallow defect states, this incor- 
rect dispersion can cause abnormally large hybridization 
with the extended band states (either conduction or va- 
lence) and subsequently partially populate these bands. 
In order to obtain accurate defect cell energies {Ed) we 
have to correct for the extra energy in the system due to 
electron populations at these higher (or lower) energies. 
For shallow donors, the correction is: 



E 



BF 



] 0(en,/c - ec){wkfn,ken,k - ec) (10) 



Where Cn^k is the k-dependent energy of state at band 
index n, ec is the conduction band minimum energy of 
the defect-free bulk after potential alignment, fn^k is the 
band occupation, and Wk is the k-point weight. O is the 
Heaviside step function. 



B. Band-gap Error 



artificially close in energy to some band states causing 
excessive hybridization and ambiguity between shallow 
or deep defect behaviour. Secondly, the range of electron 
chemical potentials used to calculate defect formation en- 
ergies will be too small, possibly incorrectly predicting 
unstable charged defect states. 

Fundamentally, the bandgap error in LDA or GGA is 
the result of a lack of continuity with respect to the num- 
ber of electrons in the exchange-correlation potentiaP^. 
This in turn leads to self-interaction error (SIE) associ- 
ated with a bias towards delocalized wave functions. The 
most common means to correct this situation is to add 
a Hubbard-like potential to penalize partial state occu- 
pancies as in GGA+U. Unfortunately, there are many 
equivalent ways to apply this correction (ie. which choice 
of 'orbitals' to apply this to) yielding the same bandgap. 
Since defect levels are sensitive not to the bandgap itself, 
but to their position relative to the host band states, this 
ambiguity can re sult in many different predictions for de- 
fect ground states^'^^. Furthermore, corrections for SIE 
are especially important for charged defect calculations 
as reducing interaction error tends to increase the ion- 
icity of the crystaP^ resulting in more ionic relaxation 
as a defect state is populated and hence greater energy 
benefit for a charged defect. 

Functionals which attempt to correct for SIE have re- 
cently emerged. One of the most robust and computa- 
tionally tractable is the HSE functionaP. The exchange- 
correlation potential is divided into short and long-range 
components via a screening length parameter. For the 
long-range portion, things are unchanged from a typical 
GGA calculation. Short-range interactions have a por- 
tion of exact Hartree-Fock exchange mixed into the ex- 
change potential, which partially corrects for SIE. This 
approach has been shown to greatly increase the accu- 
racy of the bandgap as well as relative band positions 
for a wide variety of materials^. In this work, the 
amount of short-range HF-exchange mixing is set to 25% 
which is the amount suggested by the adiabatic connec- 
tion theorem^^. We use a screening length of .1 A~^ 
instead of the more typical .2 A~^ to account for the low 
dielectric strength, and hence poor screening of ZnsP2- 
This results in close agreement between our calculated 
bandgap and experiment (see Table [l|). 

However, the added memory requirements associated 
with the use of the HSE functional for a 2x2x2 supercell of 
ZnsP2 (320 atoms) make it computationally intractable 
with our existing resources (even smaller supercells rep- 
resent a challenge). Here we have decided to use the 
'perturbation extrapolation' method put forth by Lany, 
et al.^ This is based on the idea of expressing the defect- 
influenced states in the basis of the states of the perfect 
bulk (assuming these form a complete basis): 



The most common exchange-correlation potentials, 
LDA and GGA, severely underestimates the bandgap of 
most semiconductors. This has two damaging effects on 
defect calculations. The defect-induced states may lie 



n,k 



(11) 



If we model the bandgap correction of the HSE func- 



FIG. 3. Perturbation Extrapolation Workflow - assuming that 
the perfect unit cell GGA forms a complete basis for the 
supercell defect wavef unctions, we project the GGA defect 
states onto the GGA perfect states and then use the offsets 
between the GGA and HSE unit cells to extrapolate HSE 
supercell behavior. 
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tional as a perturbation (Hp) of the perfect bulk system 
Hamiltonian (Hi^uik) via a multiplier A, the band energies 
shift: 

en,/c(A) = (^n,/c|i^6nZ/c|^n,/c) + A(^n,/c l^^p |^n,/c) (12) 



en,fe(A) = en,fe(0) + A 



5en,fc(A) 
dX 



(13) 



If we apply this same perturbation to the defect system 
(with Hamiltonian i^^nZ/c + Hd) then we have: 



FIG. 4. Tetragonal Zn^P2 has a 40-atom unit cell with 
P42/nmc symmetry. 




the tendency of the defects to not be very localized. The 
number of bands we need to calculate in the expansion 
of Eq. (11) is thus fairly limited. 

The general workflow is shown in Figure [3] Perfect 
bulk unit cell wavefunctions are determined with GGA 
and HSE functionals. Then a 2x2x2 defect supercell is 
calculated only with the GGA functional. The defect 
supercell GGA wavefunctions are then projected onto the 
unit cell GGA wavefunctions obtaining A^^u- Each band 
is then offset by the amounts prescribed by comparing the 
GGA and HSE unit cells. We can then make a prediction 
for what the supercell defect bandstructure would be if 
we had been able to utilize the HSE functional. 



III. RESULTS 



eD(A) = {^D\H,uik + Hd\^d) + \{^d\Hp\^d) (14) 



eD(A)=ei)(0) + A 



dX 



(15) 



Under the assumption of first-order perturbation the- 
ory and via Eq. (11) we have the final result: 



eD{X) = eD{0) + X^Alj 



n,k 



^e„,fc(A) 
' dX 



(16) 



As long as the assumptions of first-order perturbation 
theory are justified (ie. unchanged wavefunctions upon 
application of the band-correcting perturbation Hp) we 
would predict the defect states to track the corrections 
to the perfect bulk states in proportion to the square of 



the coefficients in the expansion {A^ ^) of Eq. (11). As 



HSE doesn't significantly affect the dispersion of bands so 
much as their relative positions to each other, we expect 
that this approximation should be justified. 

There is still a question as to whether the perfect bulk 
states form a reasonably complete basis to describe the 
defects. We are helped here by the fact that the static 
dielectric strength of Zn^P2 is low (~ 3) which leads to 



At room temperature and atmospheric pressure ZnsP2 
forms a tetragonal phase with a 40-atom unit cell pos- 
sessing symmetry P42/nmc. It is derived from the cubic 
fluorite structure with zinc at the center of a distorted 
phosphorus tetrahedra and the phosphorus surrounded 
by eight zinc sites lying roughly at the corners of a cube, 
only six of which are occupied^^ - see Figure [i] The 
calculated structural and basic thermodynamic param- 
eters are summarized in Table [H Results for the lattice 
constant and formation enthalpy are in good agreement 
with experiment with both GGA and HSE functionals. 
However, the bandgap is severely underestimated with 
GGA, which at 0.32eV is only about 1/5 of the exper- 
imental result. Using the HSE functional, we find the 



TABLE I. Calculated lattice constants, bandgap, and heat 
of formation of tetragonal Zn^P2 using both GGA-PBE and 
HSE functionals (using .1 A~^ screening length and 25% HF- 



exchange mixing). 










functional 


a (angs) 


c/a 


E, (eV) 


AHfieV) 


GGA-PBE 
HSE 


8.108 
8.160 
8.097 


1.408 
1.390 
1.286 


.32 
1.42 
1.49 


-1.79 
-1.32 
-1.53 



FIG. 5. DFT Bandstructure - calculate with HSE functional. 
Partial density of states shows makeup of VBM mainly of 
phosphorus p-character, while the CBM is mixed p and s- 
character from phosphorus and zinc states. 







FIG. 6. Slab Geometry for band alignment determination. 
Red regions show areas where the potential was integrated in 
the bulk and vacuum regions. 
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bandgap to within 5% of experiment, highlighting the 
value of using HSE for this system. 

In order to accurately describe the defect levels we 
need to determine the effect of the HSE functional on 
the bandstructure (see Figure Isl) relative to GGA. Cal- 
culating the band shifts from GGA to HSE requires more 
than simple bulk calculations, as the band energies given 
by DFT are only referenced to the average electrostatic 
potential of the simulation cell. For a periodically re- 
peated solid this is an ill-defined quantity which makes 
it impossible to directly compare band energies between 
cells with different constituents or for calculations per- 
formed with different functionals. However, if we have a 
region of equivalent potential (ie. vacuum) in both GGA 
and HSE systems we can compare the bulk potentials 
in each system to this region, and subsequently compute 
the HSE and GGA band positions relative to each other. 
Care must be taken that the solid and vacuum regions 
are large enough that the electron wavefunctions become 
negligible in the vacuum and the effects of the surface 
states are not felt in the bulk. In this work, we use a 4 
unit-cell slab along the non-polar [100] direction and an 
equal amount of vacuum (resulting in more than 35 A 
of vacuum) - see Figure [6J Non-polar, non-reconstructed 
surfaces were used to avoid creating surface dipoles with 
the accompanying undesirable step in potential across 
the solid- vacuum interface. 

The alignment of GGA and HSE bandstructures is 
shown in Figure [7| The lowering of the VBM with HSE 
is due to the reduced self-interaction of the phosphorus 
p-orbitals which primarily form the highest-energy va- 
lence bands. The upward shift of the CBM is due to 
the reduced hybridization of the Zn-s and P-s orbitals 
that make up the lower part of the conduction band with 
the valence band states. With this alignment calculated 
we can compute the GGA-to-HSE band offsets needed 
for the 'perturbation extrapolation' method discussed in 
Section II B. 

We are now in the position to apply all the corrections 
detailed in Section II. We turn our attention to the pre- 



FIG. 7. GGA to HSE Band Alignment. Bandgaps for GGA 
(left) and HSE (right) are labeled at the gamma point. Cal- 
culated offsets for the VBM and CBM are shown. 




dieted stability and doping effects of the various intrin- 
sic defects. Here it is useful to make a few notes about 
the following discussion; first, we may find that a defect 
level lies outside of the bandgap region predicting that 
this charged state is unstable as a defect localized state. 
However, it may still be possible to bind electrons to the 
defect through electrostatic interation and form hydro- 
genic effective-mass states just within the band edge. A 
discussion of these states is not in the scope of this work. 
Furthermore, we ignore the effects of both formation vol- 
ume and formation entropy in computing the defect for- 
mation energies. Formation volume is directly related to 
the change in volume when a defect is created, but is typ- 
ically only important for degenerate defect regime or for 
very high pressures. Formation entropy for point defects 
is typically on the order of a few ks so are only important 
for very high temperatures - additionally, the entropy of 
point defects in the same system tend to be similar and 
so largely cancel when we compare the likelihood of one 
defect over another. 



A. Zinc Vacancy 

Zinc is surrounded by four nearly equidistant, tetrag- 
onally coordinated phosphorus atoms. Removing a zinc 
atom from the lattice leaves four dangling bonds from the 
neighboring phosphorus atoms with mainly p-character - 
see FigurejS] These bonds are occupied with only six elec- 
trons, forming three low energy and one empty higher- 
energy defect state. This empty bonding state can cap- 
ture electrons and form an acceptor defect with a -1 or -2 
charge. As the defect states have a strong valence band 
character (e.g. in terms of symmetry), we would expect 
shallow defect behavior. 

It is instructive to study the structural relaxations 
around the Vzn site as this is intimately related to the 
occupancy and energetics of the defect state. Removing a 
zinc atom to form a neutral defect causes the four neigh- 
bouring phosphorus atoms to relax away from the defect 
site as they seek to maximize their bond overlap with the 
remaining zinc atoms in the lattice. However, this move- 
ment also shifts the dangling bonds in the vacancy region 
higher in energy, fighting this relaxation. The phospho- 
rus atom closest to the defect site (the P-atom to the 
'north' of the vacancy site in Figure Isl) has more of its 
density in the vacancy region and subsequently relaxes 
the least. As the defect becomes occupied (as in the -1 
or -2 charged states) the electrons are chiefly populating 
around the phosphorus ions which cause a further relax- 
ation away from the defect due to increased coulombic 
repulsion. As the relaxation is relatively slight in all de- 
fect charge states, we would expect that the defect itself 
would be delocalized. This is born out as the charge den- 
sity is poorly screened and well dispersed throughout the 
lattice. 

Getting the correct energetics and occupancy of the 
defect states is primarily important for their effect on 
the formation energy of the defects, which is our central 
goal. Looking at the formation energy plots in Figure |9) 
zinc vacancies exist as charged defects for all but very 
p-type regimes (fermi level close to VBM). The small 
energy difference between the different charged states 
for a fermi level at the VBM can be expected from the 
small difference in lattice relaxations associated with the 
charged states of Vzn- The formation energies in the 
n-type regime are low enough in the Zn-deficient growth 
environment for this to be an important defect. 



B. Phosphorus Vacancy 

Phosphorus is surrounded by six nearest neighbor zinc 
sites, four being almost equidistant and two zinc atoms 
about 15% further away. Phosphorus normally exists in 
a -3 oxidation state and when we remove a phosphorus 
atom the electrons from the surrounding zinc atoms have 
only high energy bonding states to move into with mainly 
Zn-p and Zn-d character - see Figure 10 We would ex- 



FIG. 8. Vzn Electronic Structure - Partial Density of States 
of the neutral defect for the four nearest neighboring phospho- 
rus atoms. Defect states highlighted in yellow, show mainly 
p-character. Bond lengths are given relative to distance to 
defect site in perfect bulk cell. 
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FIG. 9. Vzn Formation Energy. Acceptor- type defect, plotted 
for the 0,-1,-2 charged states as the fermi level varies from the 
top of the VBM (p-type) to the bottom of the CBM (n-type). 
Favorable growth condition on the left. 
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pect that these states would want to depopulate and 
form +1 and +2 charged defects as these electrons are 
no longer needed for bonding a phosphorus atom to the 
lattice. The symmetry of the defect states is not similar 
to the conduction band symmetry so we would expect 
deeper defect behavior than for Vzn- 

Removing a phosphorus atom results in a defect state 
that is effectively screened by the neighboring zinc atoms 
and highly localized between the nearest zinc sites. The 
high degree of localization results in a large lattice re- 
laxation, especially for the two zinc atoms closest to the 
defect site. Since the majority of the electron density is 
concentrated in the defect region, the four closest zinc 
atoms relax closer to each other in order to lower the en- 
ergy of the occupied defect localized states. The two fur- 
thest zinc atoms relax away from the defect region due to 
increased coulombic repulsion from the similarly charged 
zinc atoms closer to the defect. As the defect becomes 
depopulated there is a slight outward relaxation since the 
benefit in lowering the energy of the defect states is re- 



FIG. 10. Vp Electronic Structure - Partial Density of States of 
the +2 charged defect for the nearest neighboring zinc atoms. 
Defect states highlighted in yellow, show mainly Zn-p and 
Zn-d character. Bond lengths are given relative to distance 
to defect site in perfect bulk cell. 






FIG. 12. Zrii Electronic Structure - Partial Density of States 
of the -2 charged defect for the interstitial Zn site. De- 
fect states highlighted in yellow, show mainly Zn-s character. 
Bond lengths are given relative to distance to defect site in 
perfect bulk cell. 








FIG. 11. Vp Formation Energy. Donor- type defect, plotted 
for the 0,+l,+2 charged states as the fermi level varies from 
the top of the VBM (p-type) to the bottom of the CBM (n- 
type). Favorable growth condition on the left. 




duced and the zinc atoms want to increase their bonding 
with the rest of the lattice. The high degree of localiza- 
tion would also suggest deep defect behavior. 

Looking at the formation energies in Figure [TT] we 
see deep defect behavior where the defect only becomes 
charged for fermi levels in the neutral to p-type regime. 
Even for favorable growth conditions these defects are 
too high in energy to play a significant role in the 
dopability of zinc phosphide. 



atoms tetrahedrally coordinated with four phosphorus 
atoms and in-line with zinc atoms in adjacent planes. 
Incorporating an extra zinc atom into the lattice forms a 
localized defect state of mainly ionic Zn-s character since 
the neighboring phosphorus atoms have closed-shell con- 
figurations. As for zinc atoms in the perfect bulk, the 
interstitial zinc has a tendency to depopulate the s-shell 
states. We would expect donor defect behavior with a 
+1 or +2 charged state and this is what our calculations 
show. The symmetry of the defect states are similar to 
the conduction band character, consequently we would 
expect shallow behavior. 

The neighboring zinc atoms are affected the most by 
the incorporation of the Zui defect as they relax away 
from the similarly charged interstitial ion. The phospho- 
rus atoms relax closer to the interstitial site in general 
as there has been a partial charge transfer from the Zui 
defect. There is only minor differences in relaxation be- 
tween the various charged states and the defect is poorly 
screened with a delocalized wavefunction centered on the 
majority of the zinc atoms in the lattice. All of which 
would suggest shallow defect formation. 

In Figure [iSJ the formation energy plots show shallow 
behavior where the defect becomes charged for all but 
fermi levels high in the n-type regime. As the fermi level 
drops to the VBM the formation energy of the defects 
becomes very small for the favorable growth conditions 
(Zn-rich regime). These defects should be important to 
consider. 



C. Zinc Interstitial 

For both zinc and phosphorus interstitials we voxelized 
the ZnsP2 unit cell and tested the sites in order of fur- 
thest distance from neighboring atoms. Zinc interstitials 
are predicted to form in the voids on the zinc plane of 



D. Phosphorus Interstitial 

The most energetically favorable position for the phos- 
phorus interstitials are in the voids in the zinc plane of 
atoms nearly equidistant from three zinc atoms and di- 
rectly above a phosphorus atom displaced from the plane 
below. The additional electrons form states of mainly 



FIG. 13. Zrii Formation Energy. Donor- type defect, plotted 
for the 0,+l,+2 charged states as the fermi level varies from 
the top of the VBM (p-type) to the bottom of the CBM (n- 
type). Favorable growth condition on the left. 
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FIG. 14. Pi Electronic Structure - Partial Density of States of 
the neutral defect for the interstitial phosphorus site. Defect 
states highlighted in yellow, show mainly P-p character. Bond 
lengths are given relative to distance to defect site in perfect 
bulk cell. 
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covalent p-character between the Pi site and neighbor- 
ing zinc and phosphorus atoms - see Figure [T4J Due to 
the natural oxidation state of phosphorus (-3) we expect 
that the defect should be an acceptor type, which is con- 
firmed by our calculations. Since the symmetry of the 
defect states is similar to the VBM, we would predict 
shallow defect behavior. 

The neutral defect state is highly localized, so Pi causes 
relatively large ionic relaxation away from the defect site, 
though for the two zinc atoms most involved with the 
covalent bonding this relaxation is reduced because this 
movement also increases the energy of the populated de- 
fect states. The charged state is poorly screened and 
very delocalized across the lattice. Consequently, there is 
very little relaxation as the defect captures electrons from 
the conduction band. Apart from energetics, the severe 
delocalization of the charged state and good symmetry 
match between the defect and valence states causes al- 
most band-like behavior where we should readily capture 
electrons from the conduction band. 

The neutral defect has relatively low formation energy 
for P-rich growth conditions. Due to the small differ- 
ences in relaxation between the charged defect states, 
there is only a small difference in formation energy 
between them at the VBM. Thus, for fermi levels even 
modestly into the n-type regime the formation energy 
of the charged defects becomes very small making these 
critical defects to consider. 



E. CONCLUSION 

From our calculations, the likely candidate for the lack 
of n-type dopability is the Pi defects. Both of the accep- 
tor defects (Fzn and Pi) have low formation energies as 
we move the system into the n-type regime. While there 



FIG. 15. Pi Formation Energy. Acceptor- type defect, plotted 
for the 0,-1,-2 charged states as the fermi level varies from the 
top of the VBM (p-type) to the bottom of the CBM (n-type). 
Favorable growth condition on the left. 
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are more suitable Vzn sites, the Pi defect are significantly 
less costly and should be the vastly more prevalent defect. 
The Zui defects would tend to aid n-type doping as they 
are electron donors, though they have too high a forma- 
tion energy for fermi levels even moderately n-typed to 
significantly compensate for the Pi defects. 

Thus, as we dope our system with electrons we create 
a large amount of acceptor defects which act as 'electron- 
sinks' and capture mobile electrons from the conduction 
band and neutralize the doping. In fact, the Pi defect 
requires zero formation energy for fermi-levels midway 
towards the conduction band. This would 'pin' the fermi- 
level - that is, as we try to approach this level of doping a 
massive amount of Pi defects would form in the crystal. 
Since the creation of Pj defects fight n-type doping, reach- 
ing this level would not be expected. There is some hope 
in that we would anticipate the Pi defects to repel each 
other (especially the highly charged state), thus limiting 



10 



their concentration, and some n-type doping would then 
survive. Indeed, initial results show a significant penalty 
in having two fully charged Pi defects within the same 
2x2x2 supercell. However, this still corresponds to a large 
defect concentration and there is low likelihood that the 
desireable properties of Zn^P2 such as high minority car- 
rier diffusion lengths would survive in the regime where 
Pi defects saturate. 

Solving for the fermi level self-consistently using the 
above formation energy functions of the various defects, 
we predict an intrinsic fermi-level of 0.55 eV, which is 
mildly p-type. Pushing the fermi level beyond roughly 
1.05 eV becomes impossible as the number of defects ex- 
ceeds the numer of sites in the cell at this point. 

Since the Pi defects are the problematic defects, any 
means to suppress them should help the n-type doping 
issue. For extreme Zn-rich growth conditions, where the 
acceptor Pi defects are suppressed and the donor Zui are 
enhanced there may be some hope of weakly n-type ma- 



terials being formed. A suppression of interstitial defects 
in general, such as straining the crystal as it is grown 
may prove fruitful. Cluster doping^^ donor atoms with 
ones that suppress Pi formation would be another avenue 
to explore. 
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